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ABSTRACT 

The transient behavior of the phase error probability 
density function of a first-order phase-locked loop in the 
presence of additive white Gaussian noise is determined 
using numerical techniques. In addition, the modulo-2ir 
probability density function and the cycle slipping density 
function are found. The algorithm used in these analyses 
involves numerical integration of a normalized and factored 
Fokker-Planck equation. Results are shown for cases in- 
volving various signal-to-noise ratios and initial conditions. 
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I. INTRODUCTION 



A. BACKGROUND INFORMATION 

The space age spurred great interest in the development 
of phase-locking techniques since fixed tuned receivers were 
impractical for receiving signals experiencing a wide range 
of doppler shifts. Due to the inherent nonlinearity of the 
loop phase comparator, early attempts to analyze phase- 
locked loop behavior involved the use of graphical phase 
plane methods which are summarized by Viterbi [1]. The 
initial use of Fokker-Planck techniques to determine the 
steady state probability distribution of the first order 
loop phase error was accomplished by Tikhonov [2], [ 3 ]. Since 
then much work has been done using the Fokker-Planck method 
of analysis but little has been published dealing with the 
transient behavior of the phase-error process. 

Various techniques have been employed in attempting to 
statistically describe the phase-error process of first 
order phase-locked loops.- Stationary phase-error distribu- 
tions are well documented by Viterbi [1] and Tikhonov [2] 
and [ 3 ]. Transient solutions of the Fokker-Planck equation 
have been obtained by using eigenfunction expansions [4] 
and by exploiting the similarities between the Fokker- 
Planck equation and a one-dimensional heat-flow equation 
[5]. La Frieda [4] arrived at a statistical description 
of the reduced modulo-2ir phase-error process by applying 
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separation of variables. This application yielded an eigen- 
value problem which he solved numerically. Dominiak and 
Pickholtz C 5 ] arrived at a direct numerical solution to 
the transient behavior of the phase-locked loop (PLL) by 
subjecting the Fokker-Planck equation to the same numerical 
procedures as was provided for the one-dimensional heat- 
flow equation by von Neumann and Richtmyer [6], Their 
resulting partial difference equation was solved using a 
digital computer. 

This study was started in the belief that the previous 
work in the analysis of the transient behavior of the phase- 
locked loop could be expanded and improved. La Frieda’s 
technique is applicable only for the reduced modulo-2Tr 
solution and does not describe the actual behavior of the 
phase-error process. Dominiak and Pickholtz create serious 
doubt as to the completeness and validity of their technique 
in that their conlcusion of insignificant buildup of proba- 
bility mass at 4> = ±k2ir conflicts greatly with qualitative 
estimates of Viterbi [1] and Lindsey [7] and with one’s 
intuition. It was also felt that Dominiak and Pickholtz 
did not make a good choice of SNR for their numerical exam- 
ples. The stationary probability density function for their 
examples is so flat that any buildup of probability mass 
would probably be very small. 

In this paper the Fokker-Planck equation is normalized 
using several constants defined for the purpose [1] and C 5 3 - 
The resulting equation is then expressed in terms of the 
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product of two functions, one of which is indepen- 
dent of time. The other function, which is much 
smoother than the original probability density, is found 
numerically on a digital computer as a function of time. 
Multiplying the two yields the solution to the Fokker- 
Planck equation which quantitatively describes the transient 
behavior of the phase-error process of a first order phase- 
locked loop. It is found that this procedure is very effi- 
cient computationally and readily provides numerical results 
describing the phase error process, the modulo-27r phase 
error process, and the cycle slipping statistics'. 

B. REVIEW OF PHASE-LOCKED LOOPS 

A phase-locked loop is a device which automatically 
controls an oscillator or periodic function generator so as 
to operate at a constant phase angle relative to a refer- 
ence signal source. Figure 1 shows the fundamental hard- 
ware structure of the loop with the appropriate signals. 

The three major components of the PLL are: 

a) A multiplier which multiplies the input signal to 

( 

the loop and the output signal of the voltage controlled 
oscillator, i. e. , 

x(t) = 2AK^sin0 (t )cos0 ' ( t ) 

= AK^{ sin[ 0 (t )-0 1 ( t ) ] + sin[0(t) + 0 * ( t ) ] } , 

( 1 ) 
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where e(t) and e*(t) are the input signal phase and the 
voltage controlled oscillator output phase respectively. 

b) A linear low-pass filter with zero-input response 
e Q (t) and impulse response f(t) which filters out the sum- 
frequency component of x(t) while passing a filtered version 
of the difference-frequency component, i.e.. 



c) A voltage controlled oscillator (VCO), the output 
frequency of which is equal to the sum of its quiescent 
frequency with the product of its input voltage and a 
proportionality constant, i.e., 



Following through on the mathematics [1] it can be 
shown that for e (t) = 0, 



For phase error <j> ( t ) = e(t)-e'(t) and total loop gain 
K = K 1 K 2 



t 

e (t ) = e (t) + / x ( t-X ) f ( X )dX t >_ 0 . (2) 

° 0 




(3) 



d<fr» (t) = 
dt 



t 

w o + K 2 / f(t-X)AK 1 sin[0(X)-0 , (X)]dX. 



A 



d»(t) = de (t ) 
dt dt 



o 



t 

- AK / f(t-X)sin<fr(X)dX. (5) 
0 
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Letting 



then 



e 1 (t) 


= 0(t) - w t 
0 


(6) 


CD 

ro 

/“N 

ct 


= 0 ’ ( t ) - u o t. 


(7) 


<l>(t) = 


CD 

Ct 

1 

CD 

ro 

ct 


(8) 



and 



rWt'i d0 1 ( t ) t 

^ — - AK / f (t-X)sin<j>(x)dX 



(9) 



which implies the PLL model shown in Fig. 2. 

Additive noise at the input which is narrowband Gaussian 

noise n(t) with zero mean and spectral density Jjd near the 

2 

signal frequency, can be expressed as 



n(t) = JITtn, (t )sinw t + n 0 (t)coso t} (10) 

i o d o 

where n^(t) and ng(t) are independent, zero mean Gaussian 
processes of identical spectral densities, which are the 
same as the spectral density of n(t) but translated downward 
in frequency so as to be centered about zero frequency. A 
derivation similar to the derivation of (9) is made [1] with 
the noise having the form 

n'(t) = -n^(t )sin 02 (t ) + n 2 (t)cos 02 (t) (11) 

where n'(t) may be treated as a white Gaussian random 

process with spectral density ^o. This leads to the model 

2 
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of the PLL in Figure 3, [1]. The equation of operation of 
the phase-locked loop is now given by 



a<»(t ) 
dt 



d0, (t ) t 

- Kf f ( t-X ) [Asin<f> ( X ) + n'(X)]dA. 



dt 



( 12 ) 
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FIGURE 1. Fundamental Hardware Structure of the 
Phase-Locked Loop. 




dt 

FIGURE 2. Phase-Locked Loop Modeled Without Noise. 




FIGURE 3. Phase-Locked Loop Modeled With Additive Noise. 
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II. FOKKER-PLANCK EQUATION 



Let a first-order phase-locked loop have an input signal 

e i (t) of constant frequency u , i . e ., e(t)=u)t, and a voltage controlled 

oscillator with quiescent frequency u> . Additive thermal 

noise, n(t), at the front end of a receiver is a zero-mean 

Gaussian process whose spectral density is essentially flat 

over the frequency range of the receiver and may be assumed 

to be white with value ^jo. The equation describing the 

2 

loop operation [1] is from (12) 



d<ft (t ) 
dt 



(u-u> ) - K[Asin4>(t) + n’(t)] 



(13) 



where <j>(t) is the difference between the incoming signal 
phase and the VCO output phase and K is the total loop 
gain. 

It can be shown [1] that 4>(t) is a continuous Markov process. 
The noise n’(t) is a stationary Gaussian process with zero 
mean and it is nearly white over the frequency range of 
interest. Therefore, 

N 

E[n* (t 1 )n’ (t 2 ) ] = -£• 6(t 2 -t 1 ). (15) 



Under these conditions, the probability density function 
p(<fi,t) of the Markov process phase error <J>(t) must satisfy 
the Fokker-Planck equation [1], i.e.. 
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ft“ ^ = - h [ (w-w o -AKsin<}) ) p ( 4> ,t ) ] 

+ A o » 2 pC.O 

~ lr_ H 2 

with initial condition 

p(4>,0) = 5 ( 4> — Q ) 

and normalizing condition 

oo 

/ pU,t)d$ = 1 . 

— oo 

It is convenient to express (16) in normalized form. Let 



■ L = AK/4 


loop-noise bandwidth 


(19) 


a = A 2 /N B t 
o L 


signal to noise ratio 


(20) 


Y = (w-w )/^B 


L frequency ratio 


(21) 


x = HB L t 


normalized time. 


(22) 



Substituting (19) - (22) into ( 16 ) yields a normalized 
Fokker-Planck equation which is 

= fr- [(sin«j,-Y)p(*,x)] + - a2 g ( * . ^ - (23) 

3t 94> a ^2 

An analytic solution to (23) has never been found. 

It is apparent that by considering only very small phase 
errors, i.e., sin<4> - <j» , then (13) can be linearized. Making 



( 16 ) 



(17) 



( 18 ) 










* 




the same approximation in (23) gives a closed form solution 
which satisfies the initial and boundary conditions on p(4>,t). 
This solution is 




1 1/2 



exp 




which is Gaussian with 



mean = m(x) = y + (<f> - y) e 



(25) 



and 



variance = o^(t) = (1 - e ^ T )/a . 



( 26 ) 



It is the Gaussian nature of (24) which provides the 
information necessary to choose the initial boundary condi- 
tion function for the numerical algorithm used in the solu- 
tion of (23). However, before the selection of an initial 
condition is discussed, it is appropriate to analyze the 
technique utilized in arriving at the numerical algorithm 
used in this study. 
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III. TRANSIENT ANALYSIS 



A. FUNDAMENTAL TECHNIQUE 

The fundamental technique employed in this study is the 
treatment of the normalized solution to the Fokker-Planck 
loop equation as the product of two functions, i.e., 

pU,t) = gU,T) exp[a (cos<}> + y<}> ) ] . (27) 

This exponential factor was selected on both mathematical 
and intuitive grounds. Intuitively, one would expect that 
after having slipped a cycle the loop dynamics would again 
exhibit some measure of steady state behavior prior to 
another slip. Lindsey [7] makes a similar observation 
in his analysis of cycle-slipping probabilities. At any 
one time the loop phase-error tends to wobble about one of 
an infinite set of points spaced by 2 tt radians. It seems 
reasonable that the phase error probability density should 
resemble a mod-2TT stationary density function in each 2 tt 
interval with a broad taper away from the origin reflecting 
the random walk, i.e. cycle slipping, nature of the problem. 

The "envelope” encompassing this succession of quasi- 
mod-2ir functions can be easily imagined to be much broader 
and smoother than the p ( 4> , x ) function. It seems plausible 
that the function g(<j>,t) from (27) is this smooth "envelope". 
The desirability of this observation is that the "smoothness" 
of g(<t>,x) would allow a much coarser grid to be used in a 
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numerical integration scheme than would be possible with 
the function p(<}>,t) . 



Applying (27) to the normalized F-P equation (23) one 
obtains 



It is now apparent that solving for g(<J>,r) is computationally 
simpler than solving for p(<t>,-r). And, since the exponential 
factor is a constant for each value of 4> along the <f>-axis, 
p(4>,t) can be reconstructed readily using (27). 

Equation (28) is evaluated numerically by using the 
central-difference method to approximate the derivatives, 
such that with A 4> and At being the increments in <(> and t 
respectively, then 




( 28 ) 



9g(<fr g(4> , t + At )-g( 4> , t ) 



(29) 



At 



3g(<> , l ) „ g(<fr + A» , t )-g(4> , t ) 



(30) 




„ g(<j) + A<}> , T ) - 2g ( (ft , t ) + g ( <f> — A<» , X ) 
(A<fr) 2 



(3D 



and it follows directly that 



(32) 

g( 4> ,t+At )=g( ,t ) + ( Y-sin<{> )^-[g( 4 >+A<j> ,t )-g( <f> — A<f> ,t ) ] 



+ “• p [g( <t> + A(J> ,T )-2g(<t ,T )+g( <t>-A<*> ,T ) ] . 

a ( A4> ) 2 



17 



B. OTHER CONSIDERATIONS 



1. Initial Condition 

It is the response of (23) to an input impulse 
function (17) which is sought. However, it was felt that 
a finite difference method algorithm acting on the extreme 
discontinuities of a numerical delta function of A <t> width and 
1/A* height might result in significant errors in the re- 
sults. Therefore the Gaussian nature of (24) was utilized 
in choosing an initial boundary condition for the algorithm. 
The boundary condition selected was a Gaussian function 
which corresponds very nearly to the response of (23) at 
some very small time, assuming the impulse function (17) 
existed at time x = 0. To check the accuracy of this boun- 
dary condition three techniques were employed. The first 
was to plot the variance of the calculated phase error 
versus time to determine whether or not the loop is opera- 
ting in its linear region; if so, the variance vs. time 
curve would extrapolate linearly back to an equivalent 
origin, where the variance is zero at x = 0. The second 
was to actually test the phase error response to the numer- 
ical delta function mentioned above. The third was to 
compare stationary and transient mod-2Tr results with those 
already published [1], [4] and [5]. 

When the initial condition is the Gaussian function, the 
equivalent time, x , at which this condition would occur 
if p(*,0) = 6 C 4>— <J> ) is given by (26) as 
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( 33 ) 



2 2 , , _ 1 - e 2t o 

°o = ° (T o } 5 



and it follows that 



t q = - i[ln(l-aa Q 2 ) ] . 



(3*0 



In this study B^t = 0 and B^t = -jj- are the times at which 
the equivalent delta function and the Gaussian initial 
condition existed respectively. In [5], however, B^t = 0 
was taken to be the time of existence of their Gaussian 
initial condition, and this does not truly represent the 
dependency of p(<J>,t) on time. 

2. End-value Algorithm 

The use of a digital computer in evaluating this 
numerical integration scheme requires a finite interval 
which comprises the sample points along the <j>-axis. A 
potentially large source of error in any computations per- 
formed on this fixed interval is the method by which the 
end points of the interval are evaluated after each At step. 
The central-difference method algorithm in this study, if 
applied exclusively, would result in an error propagating 
inward one a<$> increment per at increment starting with the 
first at increment. This error would make the procedure 
useless for large t unless the total 4>-interval were made 
prohibitively large. To minimize error from this source. 
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a sixth-order Stirling extrapolation formula is used to 
evaluate both of the end points. 

3 . Stability 

Stability in a numerical algorithm implies that the 
errors generated at each step of the numerical integration 
process tend to decay as the solution progresses in time. 

The method used in this study to determine the stability 
conditions is patterned directly after von Neumann's Fourier 
series technique [6] using the format presented in [8], The 
mesh is assumed sufficiently small such that the coeffi- 
cients in (32) may be treated as constants over a region 
which is small but contains many mesh points. For N sample 
points spaced A<J> radians along the ^-interval, the initial 
errors, i.e. t = 0, which are assumed to be located at the 
sample points can be expressed for each point in complex 
exponential form as 

N-l jB k(A<t>) 

E, = E A e n , k - 0,1,2,. N-l (35) 

K n=0 n 



where 



n 



n2ir 

(N-1)(A$) 



(36) 



Because the finite-difference equations will always be 

linear and, therefore, separate solutions additive, one need 

only consider the propagation of the error due to a single 

term, e.g. e J&k(A<}>), ip^ e coefficient A is a constant and 

n 
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can be neglected; and for simplicity the subscript n on 
6 is dropped. 

To investigate the propagation of this error for 
increasing t, one must find a solution of the finite-diff- 
erence equation which reduces to e «JBk(A<|>) when t = £(At) =0. 

t h 

is defined as the error at the k point along the <j> 
t h 

interval at the iteration along t. Assume that the 

error has the form 



jBk(A4) Xt = j 61c ( A ) X^(Ax) 
kit e e 

= 6k( A<j) ) ^ 



(37) 



where £ = e^ Ax ^ and X, in general, is a complex constant. 
Clearly, (37) reduces to e 3Bk(A<|>) w hen l = 0. Also, the 
error will not increase for increasing t provided 



Ul i 1. 



(38) 



Since the error r, . satisfies the same difference 

k2. 

equations as g(<t>,x), substitution of (37) into (32) and 
subsequent division by e <JBk(A<|>) £ A yields 



£ 



1 + ( y-sin^ )|^-(exp[ j 6 ( A 4> ) ]-exp[-j B ( A 4> ) ] > 

+ - — ^ {exp[j 6 ( A 4> ) ] — 2 + exp[-j B ( A<}> ) ] } 
° (A*) 2 



1+2 (y-sin<}> )^-sin[B ( A<£ ) ] + (f) — 

A<fl a ( A <}> ) 2 



{cos [ B ( A<j> ) ]-l } 
(39) 



21 



If one assumes that the 2 ( y-simj) )^-sin[ 6 ( A4> ) ] term is 

A <J> 



negligible for large N and very small. A<j>, then from ( 38 ) 
the following inequality is established: 



For numerical calculations let At and A<j> be related as 



For values of a in the vicinity of unity and for A<j> << 1, 
then At << A<)>. Hence, the dropping of the term containing 



4 . Cycle Slipping 

The inclusion of absorbing boundaries around some 
sub-interval of <j>, e.g. [-ir,ir], enables one to determine the 
probability that the phase error has not exceeded the limits 
of that sub-interval at any point in time. As the density 
function spreads out, any probability mass which reaches 
the limits of the sub-interval is removed from the system 
entirely by forcing the density function to be zero outside 
the sub-interval. The effect of this technique is to make 
the remaining area of the density function equal to the 
probability that the phase error has not reached the limits 
of that sub-interval. 

The cycle-slipping density function is denoted 
q(<J>,T) to distinguish it from the corresponding function 




(*» 0 ) 



At - jj-( A <j> ) ^ . 



('ll) 



— is justified. 
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p(4>,t) for the unbounded case. As long as the phase error of 
the loop remains within the limits of a specified subinterval, 
the density function qO|>,x), with initial condition 

q(<i>,0) = 6 ( <(> ) , will satisfy the Pokker-Planck equation (23), 

i. e . 



Following the development from Viterbi [1], for y = <t> Q = 0 
the probability that the phase error 4> has not reached some 
value <t>^ by time t is given by 



The fundamental difference between q(<j>,x) and p (<|>,x) is 
that q(<j>,x) is not strictly a probability density function, 
unless it is normalized by ^ ( x ) . The numerical algorithm 
used to generate q(<f>,x) is the same as for generating p ( 4> , x ) 
except that the end points q(<f ,x) and q(-<j>^,x) are forced 
to be zero. The end points for p ( 4> , x ) are found by extrapo- 
lation, thus simulating an infinite interval. 

In this study an analysis was made of the cycle 
slipping probabilities for the primary cycle slip, i.e., 




<Kt) = q ( <j> , x )d<J> 

t 0 



(43) 



where 



qU £} x) = q(-<|> ,x) = 0 for all x . 



(44) 



= [“XT , T7 1 



(45) 
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and for the primary half-cycle slip, i.e.. 



C — 4> ^ , 4> £ 1 = C-ir/2,ir/2] . (46) 

Plots of H'(t) versus t show the time dependent behavior of 
this cycle slipping probability. 



IV. RESULTS 



The entire numerical analysis of this study was pro- 
grammed for an IBM-360/67 computer using Fortran IV language, 
level G. Double precision computation was used exclusively 
to minimize round-off and truncation errors. The initial 
condition on p(<*>,t) in (23) was a Gaussian function with 
mean <»> and standard deviation a = n/20. The ^-interval of 
investigation was [-117 t,11tt] with g(<f>,x) and, therefore, 
p(<j>,-r) being computed over the entire interval to allow for 
the non-symmetrical conditions y / 0 and <)) / 0. Stability 

considerations specify that At < ( a/2 ) ( A ) , and to accom- 
modate this requirement At was chosen such that 
2 

At = ( a/4 ) ( A 4> ) . The sampling interval along the <j>-axis 
was selected to be A 4> = tt/ 50 except in the cycle slipping 
calculation where A<f> = u/100. For the examples which will 
be shown here, it was found that negligible improvements 
were obtained for smaller values of A<p and a . In fact, 
a "delta function" of width A 4> and height 1/ A was used as 
the initial condition on p(4>,t) in several cases, and the 
results obtained varied only by approximately 1 part in 
20,000 from those obtained with the initial Gaussian condi- 
tion above. Table I indicates the various conditions which 
were studied. The specific values were chosen so as to 
provide comparison with previously published work [1], [4], 
[5]. Also included in the table are the accuracy 
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I 









considerations for this study expressed as areas of the 
probability density function after the indicated number of 
time steps. 



TABLE I 



a 


<j> 


Y 


At 


Number of 
Time Steps 


Area Under 
Curves 


iH 

• 

O 


0 


0 


TT^XlO"^ 


20,000 


1.000 


1.0 


0 


0 


Aio - 11 


40,000 


1.000 


0 

• 

1 I 


0 


sin( ir/4 ) 


2 

TT XlO 


20,000 


1.000 


2.0 


tt/2 


0 


2 -4 

2it xlO 


20,000 


1.000 



Figures 4 , 5 , 6, and 7 are graphs relating to the condi- 
tions a = 1.0 and y = <j) = 0. Figure 4 is a sample of the 
dynamics of the g(4>,x) function from (27). It demonstrates 
the "smoothness" of the function with respect to the p(4>,x) 
function which is shown in Figure 5 • Transitional and 
stationary modulo-27r densities are depicted in Figure 6. 

The stationary modulo-2iT density function is seen to be a 
reproduction of the curve determined by Viterbi [1], Of 
particular interest from Figures 5 and 6 is that the p((j>,x) 
function continues to spread out with buildups of probability 
mass at the ±k2ir points even though the steady state modulo-2ir 
density function has been reached. Figure 7 is a plot of the 
variance of the phase error probability density function 
versus time. It is interesting to note the linear nature of 
the curve for B. t > 1.2. 
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La Frieda's eigenfunction analysis [4] resulted in 
transitional and steady state mod-2ir densities which were 
reproduced using the algorithm of this study. Figure 8 
is included to show the p(4>,t) curves corresponding to 
each of the times used in [4] with an additional curve 
to demonstrate clearly the asymmetric probability mass 
buildup associated with the conditions a = 2, y = 0, and 
<j> o = tt/2. 

Inasmuch as La Frieda used a different normalizing 
scheme. Dr, in his study, equivalent times, i.e. B^t=DT/2, 
are shown on 'the graph for ease of reference. The 
asymmetry of probability mass distribution is due to the 
main body of p(<>,t) being offset from <J> = 0 for some finite 
period of time while it is shifting. Since p ( , t ) is 
settling even while it is shifting, more mass will accumulate 
at = 2 tt than at <t> = — 2 -rr . Once the main body of p(<j>,x) is 
settled about <j> = 0, however, the slippage of mass in 
either direction becomes equally likely and the level of 
asymmetry will remain constant for increasing time. 

Viterbi dealt with a detuned loop which was also ana- 
lyzed in this study. For this case a = 1.0, y = sin(7r/4), 
and <t> Q = 0 [1], Figures 9 and 10 show the mod-2-[r prob- 
ability density function and p(4>,t) for these conditions 
respectively. As before the steady state mod-2Ti density 
function determined by this study is a reproduction of that 
found by Viterbi. As would be expected and is demonstrated 
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in these two curves, the detuned loop exhibits a definite 
tendency to slip in the direction of the steady-state phase 
shift caused by detuning. 

Figure 11 is the p ( 4 > , t ) function corresponding to 
a = 0.1 and y = <$> =0. This set of conditions is one of 

those selected by Dominiak and Pickholtz [5] in their 
analysis of the first order PLL by a numerical technique. 
This curve shows the buildup of probability mass which 
went undetected in the other study. It is felt that Domin- 
iak and Pickholtz did not find this buildup for several 
reasons. Firstly, they did not look look enough in their 
time frame. Secondly, their attempt to deal directly with 
the p(4>,t) function rather than an intermediate function, 
e.g. g ( 4> , t ) in this study, resulted in a grid mesh size 
which, although mathematically correct, was too coarse to 
detect the subtle and yet significant mass slippage of the 
probability density function. Finally, from [1] it is 
apparent that for a = 0.1 the mod-27r density function is 
very flat and is not a good choice of SNR with which to 
search for probability mass buildup. 

Figures 12 and 13 show the dynamics of the cycle slip- 
ping density functions, q ( 4» , t ) , for primary cycle slips, 
[-ir,ir], and primary half-cycle slips, [ — tt/2 , tt/ 2 ] . These 
curves are plotted for the conditions a = 0.5 and $ o = y =0. 
In Figures 14 and 15 are the curves which demonstrate how 
the cycle slipping probability, ^(t), varies with time for 
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the cases shown in Figures 13 and 14. The function q(4>,-r) 
for the primary cycle slip is essentially identical to the 
p(<{>,t) function until a significant amount of probability 
mass is lost at the absorbing boundaries at <j> = ±tt. Until 
this loss occurs the area under the q(<j>,x) curve, i.e. 4>(t), 
is 1.0. The time at which ij>(t) drops below some value, 
e.g. 0.5, can be determined from Figure 14. From Figures 
13 and 15 similar data may be obtained for the primary 
half-cycle slip which has absorbing boundaries at <J> = ± ^ . 
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FIGURE g($,T) Density Function for <x=l. 0, y=$ o =0. 
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Phase Error, <j> (6.0 rad/inch) 
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(qouf/s^Tun 2*0) Joaas osaqd jo jpd - ( x c ^ ) d 
FIGURE 5- Probability Density Function of Phase Error, 

a = 1.0, Y = 4> o =0. 
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Phase Error, 4> (6.0 rad/inch) 




•FIGURE 6. Modulo-27r Probability Density Function of 
Phase Error, a=1.0,Y = <j> =0. 
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phase error, <fr (0.8 rad/inch) 
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FIGURE 
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Probability Density Function Variance as a 
Function of Time, a=1.0. 
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FIGURE 8. Probability Density Function of Phase Error, 
ct=2.0,Y = 0, <J> =tt/2. 
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phase error, <)> (3.0 rad/inch) 
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FIGURE 9. Modulo-2^ Probability Density Function of 
Phase Error, a = 1.0, y = sin(7r/4), o = 0 . 
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FIGURE 10. Probability Density Function of Phase Error 
a=l . 0 , Y = sin(7r/4), 4> o =0. 
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Phase Error, <)> (6.0 rad/inch) 
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FIGURE 11. Probability Density Function of Phase Error 

a=0 . 1 , y=6 =0 • 
j o 
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Phase Error, <)> (6.0 rad/inch) 




FIGURE 12. Primary Cycle Slip Density Function, a=0.5. 
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Phase Error, <j> (0.8 rad/inch) 
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FIGURE 13. Primary Half-Cycle Slip Density Function, a=0.5» 
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FIGURE 14. Primary Cycle Slipping Probability as a 
Function of Time, a=0.5* 
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B T t, Time (0.5 units/inch) 




FIGURE 15. Primary Half-Cycle Slipping Probability as a 
Function of Time, a=0.5> 
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3 T t, Time (0.1 units/inch) 



V. CONCLUSIONS 



This study has resulted in the development of a valid 
technique for numerically finding transient solutions for 
the Fokker-Planck equation as applied to the first-order 
PLL. That the technique works is demonstrated by its repro- 
duction of existing published stationary mod-2Tr solutions 
in [1] and transient mod-27r solutions in [4], Results 
not previously available for the phase-error probability 
density function and the cycle slipping probability are 
presented here. Most important, however, is that the 
qualitative notion of the buildup of probability mass 
around points spaced 2ir radians has been verified. 
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